Ruderman-Kittel-Kasuya-Yosida interactions on a bipartite lattice 
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Carrier-mediated exchange coupling, known as Ruderman-Kittel-Kasuya-Yosida (RKKY) inter- 
action, plays a fundamental role in itinerant ferromagnetism and has great application potentials in 
spintronics. A recent theorem based on the imaginary-time method shows that the oscillatory RKKY 
interaction becomes commensurate on bipartite lattice and predicts that the effective exchange cou- 
pling is always ferromagnetic for the same sublattice but antiferromagnetic for opposite sublattices. 
We revisit this important problem by real- and imaginary-time methods and find the theorem misses 
important contributions from zero modes. To illustrate the importance of zero modes, we study the 
spin susceptibility in graphene nanoribbons numerically. The effective exchange coupling is largest 
on the edges but does not follow the predictions from the theorem. 

PACS numbers: 75.75.+a, 71.10.Fd, 75.10.-b 
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Exchange coupling between localized magnetic mo- 
ments leads to various magnetic phases that are impor- 
tant for spintronics applications*!* 2 -' 3 - Unlike the direct 
exchange coupling, the Ruderman-Kittel-Kasuya-Yosida 
(RKKY) interaction is mediated by itinerant carriers in 
the host material and is sensitive to the low-energy prop- 
erties near the Fermi surfaced On a bipartite lattice, the 
Fermi surface is naturally nested. Thus, it is interest- 
ing to explore the interplay between the lattice structure 
and the oscillations in the RKKY interaction. A recent 
theorem^ predicts that the oscillatory RKKY interaction 
is commensurate on a bipartite lattice, so that the effec- 
tive exchange coupling between two localized moments is 
always ferromagnetic when they are on the same sublat- 
tice, but antiferromagnetic for opposite sublattices. 

Graphene-based material o 6 ' 7 ' 8 ' are ideal candidates 
for exploring these peculiar features. Firstly, unlike 
most conventional metals, the Fermi surface of graphene 
shrinks to two Dirac points with relativistic dispersion^ 
and the electronic structure is well described by the 
tight-binding model on the honeycomb lattice. In addi- 
tion, since the successful fabrication of planar graphene 
it has become more than just an academic curiosity, as 
there is great application potential for electronic trans- 
port at the nanoscal o 11 ! 12 . Of particular interest are 
studies and proposals which demonstrate applications for 
graphene nanoribbons in the field of spintronics i 13 i 14 For 
example, graphene nanoribbons with zigzag edges sup- 
port localized states at the Fermi energy and cause the 
accumulation of spin and charge polarization near the 
edges , 15 - 16 - 17 - 18 - 19 - 20 

Here we aim to better understand the carrier-mediated 
exchange coupling on bipartite lattice and use a graphene 
nanoribbon as a demonstrating example. We calculate 
the carrier-mediated exchange coupling using the real- 
time formalism and find the results do not agree with 
the theorem in Ref. [5| obtained from the imaginary-time 
approach. The discrepancy between these two methods 
is analyzed in detail. The key lies in the zero modes 
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FIG. 1: (Color Online) A zigzag graphene nanoribbon with 
width L y — 12, length L x and sublattices represented by 
squares and circles. Each site is defined by integral values 
of x and y, e.g., the rectangle contains all sites with x = 2. 



of the system which require extra care to be taken with 
the analytic continuation connecting the imaginary-time 
Green's function to the real one. To illustrate the dis- 
crepancy between the two approaches explicitly, we nu- 
merically compute the RKKY interaction on the zigzag 
graphene nanoribbon shown in Fig. [T] Our results clearly 
show that the effective exchange coupling is not always 
ferromagnetic for two moments on the same sublattice, 
and neither is it always antiferromagnetic on opposite 
sublattices. 

We start with the general hopping Hamiltonian on a 
bipartite lattice, 



(i) 



where tji = t*j to ensure hermicity, and the spin indices 
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are suppressed. For simplicity, we will only consider the 
real hopping amplitude fcy. The bipartite condition re- 
quires that tij is non- vanishing only for i,j on different 
sublattices A, B. The chemical potential [i is set to zero 
so that the system is at half-filling. Undoped graphene is 
at half-filling, and can be described by an effective tight- 
binding model on a bipartite lattice, such as that given 
in Eq. Q. 

The carrier-mediated exchange coupling is propor- 
tional to the static spin susceptibility, 



Ji 



-J 2 X ^ = 0), 



(2) 



where J is the coupling between the impurity spins and 
the itinerant carriers on the lattice. The static spin sus- 
ceptibility can be obtained from the Fourier transform of 
the retarded spin susceptibility, 

X«(t) = ^(t)(sr(t)S+(Q) - S+(0)S-(t)). (3) 

We can compute the correlation function directly in real 
time, or we can use the imaginary-time method by ana- 
lytic continuation. 

Since the hopping amplitude is spin independent, the 
correlation functions do not depend on the spin orienta- 
tions. Furthermore, the bipartite condition gives rise to 
particle-hole symmetry and leads to the useful relation, 
(ci(f)ct(O)) = e l e J (c\(t)c J {0)), where a = 1 for i e A 
and ej = — 1 for i G B. After some algebra, the spin 
correlation function takes the simple form, 



{Sr(t)Sf(0)} = e i e j Gl(t), 



(4) 



where Gy(t) = (cJ(i)cj(O)) is the single- particle corre- 
lation function. After a Fourier transform and setting 
u> = 0, the carrier- mediated exchange coupling can be 
expressed as an integral, 



J'ij J ^i^j 



dtlm[Gl(t)}, 



(5) 



where the relation Gji(-t) = [Gij(t)]* has been used. 

The single-particle correlation function can be ex- 
panded, 



Gij(t) = (4(1)0,(0)) 



X>n(O0nO>* B V(&O, (6) 



where <p n (i) is the eigenfunction with energy f„ and rip is 
the Fermi distribution function. Since the hopping am- 
plitude is real, for a given n the eigenfunctions 4> n (i) can 
be chosen to be wholly real or wholly imaginary for all 
i. On evaluating the integral one introduces the conver- 
gence factor e~' nt and 

Jij = J 2 e l t ] y^W nm (i,j)n F (i n )n F {^ m ) 2 



= J 2 €iej W nm(i,j) 



(7) 



where the product of the eigenfunctions is W nm (i,j) = 
0* (i)0„(j)0^(i)0 m (j), which is always real. Note that 
the zero modes with £ n + £ m = do not contribute to the 
RKKY interaction and that the expression in the second 
line follows in the 77 — »• limit. 

Now we turn to the imaginary-time approach. We 
choose t > so that the time ordering is straightforward. 
The spin susceptibility in the imaginary-time formalism 



Xij{r) = -{T r Sr(r)Sf(Q) 



~eiej<?y(r), 



(8) 



where the imaginary-time Green's function is Qij (r) = 
(cJ(t)cj(O)) for r > 0. In the Fourier space, the spin 
susceptibility is 

Xij (in n ) = \^ ^ dr e^Q%(r). (9) 

If the analytic continuation if2„ — > oj + irj, is performed 
as an initial step then the RKKY interaction obtained 
from the imaginary-time method is 



1 , 

'2 J 



dre^gfJr) 



Now, if i] — * the RKKY interaction becomes 





(10) 



(11) 



As the above integral is positive definite, the sign of Jij 
only depends on the product e^ . If both spins are on the 
same sublattice, eiej = 1 and the coupling is apparently 
ferromagnetic. If the spins are on opposite sublattices, 
6i€j = — 1 and the coupling is antifcrromagnetic. 

The commensurate feature described by may seem 
reasonable. For example, the tight-binding model on the 
square lattice at half filling is bipartite with a nesting 
vector Q = (tt, it). The particle- hole excitations near the 
Fermi surface carry the same momentum and produce 
an oscillatory factor cosQ • r = (— \) x+ y in the RKKY 
interaction. It is clear that the oscillation is commen- 
surate with the underlying lattice, as predicted by the 
theorem in Ref. [H. However, the analytic continuation 
itt n — > lj + irj is established in the Lehmann decomposi- 
tion, where the integration over imaginary time must be 
carried out first and the Wick rotation is the last step of 
the calculation. When the analytic continuation is per- 
formed correctly, xfji z ) — Xij( z ) where the frequency z 
is a complex number in the upper half of the plane, which 
does not include the real axis, and one can calculate the 
RKKY interaction using cither the real or imaginary time 
method. If the analytic continuation is performed before 
the integration, as is done in Eq. (fTU|) and Eq. (fTTj) . the 
solution may be incorrect. 

Let us elaborate on the results obtained from Eq. QT 
The Greens' function takes the form, 



Gij (t) 



&{i)(l>n(j)nF(£n)e 



Tin 



(12) 
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and after integrating over the imaginary time, the ex- 
change coupling is 

J 2 ee ■ 

Jij = — ^ ^ W nm (i,j)n F (£ n )n F (£ m ) 

1 _ e P(t n +U+iv) 



£n + Cm + iV 



(13) 



The sum can be separated into two parts, the first con- 
tains contributions from zero modes with £„ + £ m = 0, 
and the second contains all other terms with £„ =/= 0. 
The first part can be shown to equal 

jO = _J^l W nm (i,j)pn F (Z n )n F (Z m ). (14) 

t n +U=o 

As for the second term, the particle-hole symmetry en- 
sures that each energy £„ is paired with an equal energy 
of opposite sign £ fi = — £„, and that their eigenfunctions 
satisfy <f>j%(i) = ^i4>n(i)- In addition, particle-hole sym- 
metry provides the useful relation n F (^ n )e^ n = n F (£n). 
Therefore, 



- ^2 W n m(i,j)n F (^n)n F (^ m 



o/3(«»+Sm) 
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17] 



(15) 



and it can be seen, after exchanging the dummy variables 
(n,fh) with (n,m), that the £ m + £ n ^ part of Eq. 
([TO]) is equivalent to Eq. ([7]). Therefore we have shown 

Jij = Jij + Jij- 

The discrepancy between Jy and J7ij always exists at 
finite temperature and, as we will show in the graphene 
nanoribbon example, there are cases where zero temper- 
ature also leads to a discrepancy. We have performed nu- 
merical calculations (not shown here) on simple bipartite 
lattices and verified that the RKKY interaction need 
not follow the sign rule obeyed by Jij. However, there 
are special limits where the correction jQ vanishes. If 
all single-particle states are gapped (£ m = — £ n ^ 0) and 
the temperature approaches zero, 

f3n F (t n )n F (-Cn) ~ pe~ pA 0. (16) 

Therefore, the difference between and Jij vanishes 
for a gapped system at zero temperature. 

In summary, while X«( w ) an d Xiji^n —> w + if)) are 
identical and either can be used to evaluate the RKKY 
interactions, one must be careful with the analytic contin- 
uation when using the imaginary time formulation. The 
analytic continuation must be taken after the imaginary 
time integration, and not before, as these two operations 
do not in general commute. If done correctly, the RKKY 
interaction will not have any contributions from the zero 
modes of the system, £„ + £ m = 0. However, taking the 
analytic continuation prior to evaluating the imaginary 
time integral will results in a zero- mode term jQ . 



While J^ will always be non-zero at finite tempera- 
tures, for systems with a regular density of states near 
the Fermi surface, the contributions from the zero modes 
scale to zero in the thermodynamic limit. One can un- 
derstand this by considering how the terms in Eq. (|14[) 
scale with volume V . The eigenf unction term scales as 
W nm ~ (1/V) 2 . On eliminating the constrained variable 
m, the summation over n scales as V J d d k in the ther- 
modynamic limit. Thus, the correction term J7fl scales 
as (1/V) 2 ■ V — 1/V and vanishes in the thermodynamic 
limit. However, there are a number of systems which do 
not have a regular density of states in the thermodynamic 
limit. Graphene nanoribbons, for example, have a singu- 
lar density of states at the Fermi surface and therefore 
JP 4 is never zero. 

To compute the RKKY interaction for the zigzag 
graphene nanoribbon, shown in Fig. [TJ we first need 
to obtain all eigenstates. Making use of the transla- 
tional invariance along the ribbon direction, we take 
the partial Fourier transform of x to k x . The opera- 
tors which diagonalize the Hamiltonian are defined by 
"0(p) = Y, y <P*(P,y) c ( k x,y) where p = (k x ,p). The solu- 
tion for \k x \ < ko = 2 cos -1 [L y /2(L y + 2)]) is 



<t>(p,V) = ±sin[py 
j,y) = sm\py], 



tp(p)], for odd y, 

for even y, (17) 



where <p(p) = tan 1 [tan(p)(i' — t)/(t' + t)] with range 
-tt/2 < <p < tt/2, and t' = It cos(fc x /2). As L y 
is finite we have the constraint (L y + l)p + tp(p) = 
TO7T with m = 1, 2, . . . Ly/2. The energy is E(p) — 
T \A' 2 +t 2 + 2t't cos(2p) . Therefore we have L y /2 posi- 
tive bands and L y /2 negative bands each associated with 
one of the L y /2 values of m. For \k x \ > ko the solutions 
for m = 1, 2, . . . (L y — 2)/2 do not change but m = L y /2 
now describes localized edge states with p = tt/2 + ry. 
Substituting this imaginary p into the \k x \ < ko solution 
with m = Ly/2 gives 



<t>{p,y) 



i-ny /2 



±e 



- „i7r(3/-l)/2 



sinh[7(L y ■ 
sinhfyy], 



y)], for odd y, 
for even y, (18) 



with (L y + 1) 7 = tanh 1 [- tanh(i)(t' + t)/(t'- t)] 
and the dispersion for the edge state is E(p) = 
-f^/t' 2 + t 2 — 2t't cosh(27). Note that the eigenfunctions 
of the edge states rapidly decay into the bulk. 

The translational invariance of the static spin suscep- 
tibility x{ x i V'i x \ y') implies that it is dependent on the 
horizontal distance \x — x'\. However, a lattice point's x 
value may not represent its true position and so we in- 
troduce the offset 5 where 8 = 1/2 for y = 1, 0mod4 and 
6 = for y = 2, 3 mod 4. Thus, the spin susceptibility 
depends on Ax + A8 = (x — x') + (5 — 8'). The numer- 
ical results at half-filling (/i = 0) and with slight doping 
(fx = 0.01) for the symmetric (L y — 12) and asymmetric 
(L y = 10) lattices are summarized in Fig. O 

Since the spin susceptibility is largest on the edges, we 
will focus on these lattice sites. When L y /2 is even, the 
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FIG. 2: (Color Online) The static spin susceptibility 
x[( x , 1); { x 'iV')] f° r a z ig za g graphene nanoribbon at zero tem- 
perature. Since it is dominated by the edge response, we 
choose y — 1 and y — L y and plot the susceptibility versus 
the horizontal distance Ax + AS for a: symmetric nanoribbon 
L y — 12 (a) at half filling fi = 0, (b) at slight doping \i = 0.01; 
asymmetric nanoribbon L y — 10 (a) at half filling /i = 0, (d) 
at slight doping /i = 0.01. 

upper and the lower edges are symmetric and the spin 
susceptibilities for y = 1 with y' = 1 and y' = L y are the 
same except for the sign. A perfect reflection symmetry 
is clearly seen in Fig. [2Ja). Note that for y = y' = 1 the 
spins are on the same sublattice (labeled as red circles 
in Fig. [T]), but the RKKY interaction = —J 2 xfj is 
not always ferromagnetic, as shown in Fig. [2a). Sim- 



ilarly, for y = 1 and y' = L y the spins are on opposite 
sublattices, but the RKKY interaction is not always an- 
tiferromagnetic. In fact, both curves show non-trivial 
oscillations which pass through zero. This provides clear 
evidence that the theorem given in Ref. |5| is incorrect. 
With only slight doping (// = 0.01), we lose the particle- 
hole symmetry and conventional RKKY oscillations ap- 
pear, as shown in Fig. HJb). 

For the asymmetric nanoribbon with L y /2 an odd in- 
teger, the reflection symmetry between the spin suscepti- 
bilities for y = 1 with y' = 1 and y' — L y is lost, as clearly 
shown in Fig. [DJc). However, the RKKY interaction for 
spins on the same sublattice (y = y' = 1) is still not con- 
sistently ferromagnetic. Upon slight doping (fi = 0.01), 
conventional RKKY oscillations again appear, as shown 
in Fig. HJd). In the asymptotic regime, the spin suscep- 
tibilities for y' = 1 and y' — L y are mirror images of each 
other due to a half lattice shift between the upper and 
the lower edges of the asymmetric nanoribbon. 

In conclusion, we show that the RKKY interaction on 
a bipartite lattice is not commensurate with the lattice, 
although it has been predicted to be so. The discrep- 
ancy arises when the analytic continuation is performed 
incorrectly, which leads to errors which can be attributed 
to the zero modes in the system. The zigzag graphene 
nanoribbon is studied numerically as an illustrating ex- 
ample which highlights the discrepancy. 
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